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^ . Abstract 

OO 

C^ ' The translational invariant formulation of the coupled-cluster method is presented 

here at the complete SUB(2) level for a system of nucleons treated as bosons. The 
correlation amplitudes are solution of a non-linear coupled system of equations. These 

r — I equations have been solved for light and medium systems, considering the central but 

^\ • still semi-realistic nucleon-nucleon S3 interaction. 
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1 Introduction 

The coupled-cluster method (CCM) was invented many years ago by Coester and Kiimmel 
P, |2l to calculate the ground-state energy of closed-shell nuclei. Since then the CCM has 
been developed for applications in a wide variety of fields, including quantum chemistry, 
quantum lattices, electron gas, or nuclear physics. Although the formalism was initially 
developed to calculate only the energy of the ground state of a closed-shell nucleus, the 
CCM is not just a theory for the ground state of closed-shell systems. Further extensions 
allow to calculate also excited states, open-shell systems, or finite temperature systems, to 
mention only a few developments of CCM (see general references in Refs. |3|, ^). Currently, 
the CCM provides a widely applicable formalism to deal with the general quantum many- 
body problem at a completely microscopic level. It works for systems of both bosons and 
fermions, regardless of the type and range of the interaction. 

The CCM wave function is described in terms of the action of a cluster operator on 
some reference function, which takes into account the required symmetry properties of 
the system under consideration. The cluster operator is sum of n-body operators, which 
create n particle-hole pairs on the reference function. In the application to finite systems, 
one faces with the well-known center-of-mass motion problem. The proper treatment of 
the center-of-mass in the framework of CCM was initiated some years ago [^, ^] at the 
so-called SUB(2) level of approximation, which is limited up to two-body operators. It 
has been shown in these references that the center-of-mass is properly removed using 
appropriate combinations of one- and two-body operators, and describing the reference 
function in terms of single-particle harmonic oscillator wave functions. Such reformulation 
of the CCM is called translationally invariant coupled-cluster (TICC) method, and it has 
been applied to the ^He nucleus, considered as a system of four bosons. For the simple 
interactions considered in Ref. |^], it was shown that the results at the SUB(2) level 
(TICC2 henceforth) may differ in a few MeV from the essentially exact diffusion Monte 
Carlo results. 

The CCM is most naturally formulated in the occupation number representation. How- 
ever, very slow convergence with respect to the cut-off on the single particle basis has been 
found. This led to consider a linearized form of the TICC2 approach, which can be easily 
converted to coordinate representation. The pair correlations are considered through a 
function depending on the relative coordinate of a single pair. The ansatz for the wave 
function of a fermionic system is 

I*) = Yl (f'^i^ij) + f'^(.rij)i^i ■ ^j) + fT{rij){Ti ■ Tj) + far(.rij){ai ■ aj){Ti ■ Tj)) |$), 
i<j 

where |$) is the reference state. The four unknown two particle correlation functions 
may be determined by minimizing the expectation value of the hamiltonian, and a single 
correlation function appears in the case of a system of bosons. This so-called transla- 
tionally invariant configuration interaction method was applied in Ref. Q to calculate 
at the SUB (2) level (or TICI2 approach) the ground-state energy of some light nuclei in 
the Ip-shell up to ^^O. The main conclusion of this work is that the TICI2 methodology 
provides a very reasonable starting point for the calculation of the binding energies of 
light-to-medium nuclei. For interactions and correlations of the V4 form, as displayed in 



the above equation, the TICK results are in suitable agreement with the ones provided 
by other methodologies. 

In order to get a better description of finite nuclei there are however some obvious 
extensions to be performed. Among them are the inclusion of all the correlations at the 
SUB(2) level, i.e. the TICC2 approximation. In Refs. |5|, ^ it has been shown that the 
additional terms give a relatively small contribution to the ground-state energy of the 
He nucleus. The purpose of the present work is to discern the importance of the non- 
linear terms beyond TICK in the description of a bosonic system. In Section 2 we briefly 
summarize the main findings of previous works [g, |^, generalizing at the same time the 
results for A^ bosons. In Section 3 the detailed equations are presented and discussed. An 
application example is given in Section 4, where a simple nuclear interaction is used to 
determine the ground-state energy of a system of nucleons treated as bosons. Finally, in 
Section 5 some concluding remarks and an outlook for further studies are given. 

2 The Translationally Invariant SUB (2) approximation for 
bosons 

The basic CCM ansatz is to write the A-body wave function as 

|*)=e^|$), (1) 

where |<I>) is a reference state, incorporating the statistical and symmetry properties of 
the system, and 

N 

S=J2Sn (2) 

n=l 

is the cluster correlation operator, sum of operators of the form 

^^ = T1^ Yl (/'i'--->Pn|5'n|0, ...,0)a+ ...a+a[f, (3) 

V'^-^ Pl,...,Pn 

whose action is to create n particle-hole excitations on the reference function. The 
SUB(2) level approximation of CCM consists of neglecting the excitations of more than 
two particle-hole pairs, so that in the initial ansatz one has S = Si + 82- 

In order to get rid of the center-of-mass problem, one could directly choose a transla- 
tionally invariant reference state, but this will spoil the meaning of the Fock space operator 
as single-particle operators, due to the center-of-mass constraint. The alternative followed 
in Ref. ||g] was to build up the reference state from harmonic-oscillator (HO) single-particle 
wave functions, which can be defined through creation operators acting on the vacuum, 

\nlm) = a^iJO), (4) 

and to write the reference state of an A-boson system as 

1^) = ^(4oo)^|0). (5) 



It is well-known that this HO uncorrelated wave function is the simple product of a Is HO 
wave function for the center-of-mass and an intrinsic wave function. Although the reference 
state is not translationally invariant, this factorization property allows one to remove the 
center-of-mass motion unambiguously. Dressing the reference state with correlations has 
the danger of spoiling the factorization property, unless the correlations are translationally 
invariant. Therefore our immediate goal is to construct excitations on |<1>) without exciting 
the center-of-mass, so that in this way it can be properly removed from the SUB(2) wave 
function. 

With only 5*1 cluster operators it is not possible to get translationally invariant corre- 
lations. However, appropriate admixtures of Si and 5*2 can do it. Let us indicate by S"*-^'^) 
the corresponding new cluster operator. Translational invariance may be imposed by: i) 
recoupling the product of single-particle harmonic-oscillator states into sums of compara- 
ble products of oscillator states for the relative and center-of-mass motion of the pair, and 
ii) imposing that the center-of-mass state of the destroyed pair of particles in the occupied 
subspace is the same as that of the created pair of particles in the unoccupied subspace. 
Since we are using HO wave functions, this is accomplished by using the Brody-Moshinsky 
brackets [H]. In the case of an N boson system, one can write the cluster operator as 

°° r nO 

5(i'2) = ^ s{n) J2 (^0, 00, 0|ni/, n^l, 0) [a+ , x a+ , ^ a^^o, (6) 

n=l ■ni,n2,l 







where (nO, 00, 0|ni/, n2/, 0) is a Brody-Moshinsky coefficient, which limits the sums to the 
cases rii + n2 + I = n. This expression reflects that two bosons in occupied states of |<l>) 
are destroyed, and the created pair of particles has zero angular momentum, its center- 
of-mass being in the HO state with all quantum numbers equal to zero. Note that the 
term with both (ni,/) = (0,0) and {n2,l) = (0,0) is excluded, as it simply reproduces 
the uncorrelated state |$). However, the terms with either {ni,l) = (0,0), n2 ^ or 
ni 7^ 0, (n2,0 = (OjO) ^^^ included. These terms give precisely the Ip-lh excitations in 
the admixture. 

It is worth noting the enormous simplification imposed by the translational invariance, 
because we have a single correlation amplitude 5(n), whose argument is just a single 
number which counts the number of oscillator quanta (divided by 2) globally excited. 
By contrast, in the unrestricted SUB(2) approximation, one deals with two correlation 
amplitudes 5i(n) and 52(ni,n2,/). 

It is convenient to use the following simplified notation for the correlation operator: 

5^^'^) = s{p)a^a^al + s{p, g)a+a+ao, (7) 

where a sum over the particle indices p, q is to be assumed. These indices refer to the 
HO quantum numbers: p = {np,lp,mp), and they are different from (0,0,0), since they 
correspond to non-occupied states. For convenience, we shall also use the notation = 
(0,0,0). The coefficients are defined as 

s{p) = 6{lp, 0)5{mp,0)2{np0, 00, 0|npO, 00, 0)S{np), (8) 

s{p,q) = S{lp,lg)6{mp,-mg)C{lp,lp,0;mp,-mp,0) 

{up + nq + IpO, 00, 0\nplp, Uglp, 0)S{np + Ug + Ip), (9) 



where C{li,l2,h','mi,m2,ms) is a Clebsh-Gordan coefficient. 

Finally, let us note that simply because S^^'"^' is invariant under translations we are by 
no means guaranteed that so does the exponential form e ' . This can be easily seen by 
considering the simple example of the squared correlation operator S'^^'^'S''-^'^'. Indeed, 
one of the destruction operators oq of one factor may be contracted with the creation 
operator aj" of the other factor, resulting in a term which is not translationally invariant. 
As was shown in Ref. [^ the unwanted terms are easily excluded by the simple device of 
taking the e"^ ' operator in normal-ordered form, so that the TICC2 ansatz for the wave 
function can be finally written as 

|*)=:e^^'''' :|$). (10) 

In principle, for a system of bosons the normal ordered form does not imply a particular 
technical problem (contrary to the fermionic case). However, the usual COM equations 
are more involved, and we found it more convenient to derive the COM equations in a 
different manner, as it will be shown in the next section. 



3 The TICC2 equations 

The usual way to solve the Schrodinger equation 

i:f|^) = £;|^) (11) 

within COM involves transforming it first to 

e-^He^\<^) = E\<^), (12) 

i.e. to define a similarity transformed hamiltonian. The motivation to use this equation is 
that the transformed hamiltonian e~'^He^ may be reduced to a set of nested commutators, 
which may result in a simplification of the algebra. Afterwards, Eq. (|l^ ) is projected onto 
Op-Oh, Ip-lh, . . . states to get algebraic equations for the amplitudes. If we were projecting 
up to A^p-A^h, the corresponding set of algebraic equations would be equivalent to the 
Schrodinger equation. However, if the cluster operator is truncated at a given SUB(n) 
level the projection up to A^p-Ah states will produce redundant equations, so the current 
practice is to project up to the np-nh states. It is worth noting that if the cluster operator 



is truncated at the SUB(n) level none of Eqs. (|T^) or (12) may be exactly satisfied, and 
moreover one should not expect in general to obtain the same results from both equations. 
The similarity transformation does not appear convenient in our formalism with the 
normal ordering prescription on exp(S'). Certainly, the inverse of the operator : exp(S') : 
is not : exp(— S) :, but it is possible to find a new operator T such that the equation 

T : e^ : 1$) = |$) (13) 

is satisfied. For the specific operator S^^'"^' given in Eq. (0), there results 



r ~ 1 - (A - l)s(p)a+ao + 



-s{p,q) + ^{N + 1){N - 2)s{p)s{q) 



a+a+al + ..., (14) 



up to two creation operators. This operator T is not the inverse of : exp(S''^'^)) : because it 
is defined when acting on the reference state |$), and this fact is reflected in its dependence 
on the number of bosons of the system. The structure of T depends on the reference state, 
because also normal ordering depends on such state. 

Another much simpler alternative to obtain the equations for the amplitudes is to 



project directly the Schrodinger equation with the ansatz (|10D 



//:e^'^•'^|cI>)=i?:e^'^•'^|$) (15) 

onto the states with Op-Oh, Ip-lh and 2p-2h. 

For hamiltonians with two-body interaction potentials like 

H = ^{oi\K\(3)d^ap + - ^ (ai,a2|^|/3i,/32)a+ a+^a/jaa/ji (16) 

one requires the expansion of : exp(5^^'^-*) : up to terms involving at most four particle 
creation operators. In terms of the quantities s{p) and s{p,q) defined in Eqs. (§) and (^, 
respectively, one may effectively do the substitution 



: e 



1 + s{p)a^a^al + s{p, g)a+a+ao + -s{p)s{q)a^a+{a^faQ 



+s{p)s{q,r)apa^a^aQaQ + -s(p)s(g)s(r)ap a^a+(a^)^ao 

s{p, q)s{r, s)a^a+a^ataQ + -s{p)s{q)s{r, s)a^a^a^at{a^fal 



+^s{p)s{q)s{r,s)a];a+a+a+{a^) a^ + . . . , (17) 

which is no more complex than the usual nested commutators. Note the reorderings of ao 
and Qq related to the normal ordering prescription. 

The projection onto Op-Oh is simply the projection onto the bra state (<1>| defined by 
the reference state, and it gives the ground-state energy: 



^=($|// :e-^*'''' : 1$). (18) 



Using the following notation: 



k{a,f3) = {nJama\K\nplpmj3), (19) 

for the matrix elements of kinetic and potential energies taken between HO single particle 
wave functions, we find 

E = iViA;(0,0) +iV2w(0, 0,0,0) 



+ N2 



k{0,p)s{p) + 2{N - l)v{0,0,p,0)s{p) + 2v{0, 0,p,q)s{p,q) 



+ N^viO,0,p,q)sip)siq), (21) 



where we have also defined 



AT. = N{N -1)...{N -i + l). 



(22) 



Again, a sum over repeated particle indices p, q is to be understood in Eq. (21). 

The ground-state energy is thus given in terms of the correlation amplitudes S{n) im- 
plicit in s{p) and s{p, q). To determine these amplitudes, we must project the Schrodinger 
equation onto lp~lh and 2p-2h states taking care that the center-of-mass state is the same 
as in the reference function. The obvious state to project on is the appropriate lp~lh and 
2p-2h admixture which maintains the factorizability of the center-of-mass in the reference 
state. Such an admixture is already known, after the discussion of Section 2: it must have 
the same structure as S^^'"^' , but without the unknown amplitudes S{n). It seems conve- 
nient to classify the excitations in terms of the number of oscillator quanta 2Nx, where 
Nx will run from 1 up to some sufficiently high value N^^x so as to reach convergence. 
In this way one obtains A'^max non-linear equations to determine the unknown amplitudes 
S{n). Once these amplitudes are known, the energy is easily computed. 

Consequently, the excited bra state required for the proper projection can be written 



as 



($1 c{m){a^) aoam + €{171,71) (a^) amar. 



(23) 



where the excitation amplitudes have a structure similar to the cluster operator, that is 
c(m) = (5(/™, 0)<5(m™, 0)2(iV^0, 00, 0|n^0, 00, 0), (24) 

c{m, n) = 5{lmJn)S{mm, -mn)C{lm, lm,0; mm, -mm, 0)(Aa;0, 00, Q\nmlm, nJm, 0). (25) 

The indices m, n refer to excited states, and cannot be identical with (0,0,0). Moreover, 
they are external indices, tied to the excitation quantum number N^, so that there is no 
summation over them. 

Thus, the equation to determine the correlation amplitudes can be written as: 



H : e^'^'^' : m 



($1 c{m){a^faoam + c{m, n){a^faman 

= E{^\ c(m)(a([)^aoam + c{m, n){aQfaman 



: e 



5(i>2) 



m. 



(26) 



Evaluating all contractions contained in this equation is not a simple task, because of 
the very large number of terms which appear after the expansion of the operators. We 
have found of great help the use of the algebraic language REDUCE|^. The core of the 
evaluation is the substitution 



axtty -^ a^ax + 5{x,y), 



(27) 



which is controlled by a let rule. This automatically reorders the full expression with 
all annihilation operators on the right. Afterwards the Kronecker delta functions are 
conveniently eliminated. 
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The left-hand side of Eq. (26) can be written in the form 



($1 c{m)(a^faQam + c{m,n){a^fama„ 
N2c{m)k{m, 0) 



H :e' 



5(i'2) . 



:|$) 



+ N2 



2{N - 2)c(m, n)k{0, 0)s{m, n) + {N - lfc{m)k{0, O)s(m) 
+2(iV- l)c{m,n)k{m,0)s{n) +2{N - l)c{m)k{0,p)s{m,p) 
+4c{m, n)k{m, p)s{n, p) + {N — l)c{m)k{m,p)s{p) 

+ A^4 (N — 2)c{'m,n)k{0,0)s(m)s{n) +2c{m,n)k{0,p)s{m,n)s{p) 
+Ac{7n, n)k{0, p)s{m, p)s{n) + {N — 1) 0(111) k{0,p)s{m)s{p) 
+2c{m, n)k{m, p)s(n)s(p) 



N2 
N2 



c{m, n)k{0,p)s{m)s{n)s{p) 

2{N - l)c{m)v{m, 0, 0, 0) + 2c(m, n)v{m, n, 0, 0) 

(A^ - lf{N - 2)c{m)v{0, 0, 0, 0)s(m) 



+2(iV -2){N - 3)c(m, n)v{0, 0, 0, 0)s(m, n) 

+4(iV - l)(iV - 2)c(m, n)t;(m, 0, 0, 0)s(n) 

+4(iV - l)(iV - 2)c{m)v{0,0,p,0)s{m,p) 

+4(A^ — 1) c{m)v{'m, 0,p, 0)s{p) + 16{N — 2)c{m, n)v{n, 0,p, 0)s{m,p) 

+4(A^ — l)c{'m)v{m, 0,p, q)s{p, q) + 4(A^ — l)c{m, n)v{m, n,p, 0)s{p) 

+4c(m, n)v{m, n, p, q)s{p, q) 

+ N4 {N -2){N -3)c{m,n)v{0,0,0,0)s{m)s{n) 

+4(A^ — 3)c(m, n)f (0, 0,p, 0)s{m, n)s{jp) 

+8(A^ — 3)c(?Ti, n)f (0, 0,p, G)s{m^p)s{n) 

+2{N -1){N - 2)c(m)v{Q, 0,p, 0)s{m)s{p) 

+4c(?TT., n)f (0, 0,p, q)s{m, n)s{p, q) + 8c(r?T,, n)v{Q, 0,p, q)s{m,p)s{n, q) 

+8{N — 2)c{m, n)v{n, 0,p, Q)s{m)s{p) + A{N — l)c{rn)v{Q, 0,p, q)s{m, q)s{p) 

+2{N — l)c(m)t;(0, 0, p, q)s{m)s{p, q) + 2{N — l)c{m)v{7n, 0, p, q)s{p)s{q) 

+ 16c{m, n)v{n, 0,p, q)s{m, q)s{p) + 8c(m, n)v{n, 0,p, q)s{m)s{p, q) 

+2c(m, n)v{m, n, p, q)s{p)s{q) 

+ A^6 2{N - 3)c{m,n)v{0,0,p,0)s{m)s{n)s{p) 

+ {N — l)c{m)v{0, 0,p, q)s{7n)s{p)s{q) 

+8c(m, n)f (0, 0,|), q)s{m,p)s{n)s{q) + 2c(?7i, n)w(0, 0,p, q)s{m, n)s{p)s{q) 



+2c(m,, n)v{0, 0,p, q)s{m)s{n)s{p, q) + Ac{m, n)v{n, 0,p, q)s{m)s{p)s{q) 



+ Ns 



c{m, n)v{0, 0,p, q)s{Tn)s{n)s{p)s{q) 



(28) 



In the right-hand side of Eq. (p^) one has 



($1 c{m){aQfaoam + c{m, n){aQ)'^aman 



: e-^^'^' : |$) 



= N2 

+ iV4 



(A^ — l)c{m)s{m) + 2c{m, n)s{m, n) 
c{m,n)s{m)s{n) 



(29) 



Let us recah that in these equations any index m, n, p, or q must be interpreted as a set 
of three quantum numbers n = (n^, ln,Tnn) of the HO states, none of them being (0,0,0). 
Let us insist again that in these expressions, indices p, q refer to particle states and a sum 
on them must be understood, indices m, n are external, and are related to the quantum 
number Nx characterizing the excitation. 

We have stressed the factors A'^2, -^4, -/Ve, and N% accompanying the correlation of one, 
two, three, and four particle-hole pairs. The results of Ref. [|| showed that in the case 
of ^He the quadratic terms give a relatively unimportant contribution to the ground-state 
energy, but we cannot assume this to be the case for larger systems. Firstly, it is clear 
that only two pairs can be created for A^ = 4 particles, and consequently the TICC2 
approximation for ^He contains only linear and quadratic terms. But more pairs can be 
excited for A^ > 4, and we see that up to four pairs are present when the interaction 
is a two-body one. Secondly, the expressions given above for the projection show that 
these terms are multiplied by increasing powers of N: quadratic terms are multiplied 
by a factor A^4, cubic terms by Nq, and quartic terms by N^. Although the amplitudes 
S{n) may depend on N one cannot discard the possibility that the non-linear terms could 
become relatively more important for high values of A'^. 

We have to solve a non-linear (actually quartic) system of coupled equations to deter- 
mine the amplitudes S{n), which can be cast in the following schematic form: 



J * max 

Y^ F{Nx,ii,i2,i3,H)<S{ii)S{i2)S{i3)S{i4) 

Hii2j«3i*4=0 



0, 



(30) 



with the convention 5(0) = 1. The number A'^j; = l,...,A'max refers to the excitation 
considered. The function F{Nx,ii,i2,i3,i4) is obtained from the expressions given above 
for the projections onto the excitations up to 2A'3; quanta. The indices i in the correlation 
amplitudes refer to the number 2i of HO quanta characterizing the excitation, which we 
limit to a maximum number A^max- The independent term is given by F{Nx,0, 0, 0, 0). The 
linear version of the present development, or TICI2 approximation, amounts to consider in 
Eq. ( |30| ) only the function F{Nx, ii,0, 0, 0) (including ii = 0), and the problem is reduced 
to solve a generalized eigenvalue equation. 

There is still a long way from Eqs. (P^) - (|29| ) to the practical equations. Rotational 
invariance of the hamiltonian merged with the specific form of s{p) and s{p, q) results in a 



drastic reduction on the number of summed up indices. Additionally, the potential matrix 
elements are conveniently transformed into Talmi integrals, referring only to the relative 
coordinate of the pair. The final formula is too long to be included here. A part of the 
non-linear equations, up to quadratic terms in the amplitudes, may be found in Appendix 
A of Ref. II, for the A^ = 4 case. 

4 An application example 

We consider a model system consisting of N nucleons treated as bosons interacting via the 
Wigner part of a nuclear interaction. These systems collapse for large values of N, and so 
one could reasonably expect that the effects related to the correlations will be enhanced. 
We have used the Wigner part of the Afnan-Tang nucleon-nucleon potential S3 [^ , which 
is considered as a semi-realistic interaction having a quite strong core at the origin. This 
will enhance again the role of the correlations. 

The price to be paid when going from the configuration interaction TICK scheme to 
the full TICC2 one is that instead of solving a generalized eigenvalue problem as was 
done in Ref. ||^ , one has to solve the non- linear system of equations (^) coupled to the 
equation determining the energy E = {^\H : exp(S'(^'^') : |$). The obvious steps are 
first to obtain explicitly the matrix of coefficients F{Nx,ii,i2,i3,i4) and afterwards to 
solve for the amplitudes and the energy. The equations depend on the harmonic oscillator 
parameter a = {muj/fiY''^ used to define the Fock space, so it is convenient to adjust this 
Tiuj parameter like in a variational calculation. 

As we have already mentioned, we may limit the number of equations by means of the 
index A^max, which measures the number of excitation quanta (2Amax^w) related to a given 
term of the correlation. One may then proceed sequentially, starting from a small value of 
Arnax, and taking its solution as the approximate solution corresponding to the next set 
of Amax + 1 equations. The Newton-Raphson method turns out to be very efficient once 
we are quite close to the solution. 

An alternative way is to introduce in Eq. (30) a quenching factor q as follows: 
^ F(A,,n,0,0,0)5(ii) 

Ji=0 

-^'inax -'''max 

+ *? E E F{N^,ii,i2,iz,H)S{ii)S{i2)S{h)S{H)=^, (31) 

«i=0 (i2,i3,M)^(0,0,0) 

where in the second sum the indices (i2 5^3;^4) cannot be simultaneously equal to zero, as 
this case is included in the first sum. We let the quenching factor q slowly grow from to 
1, in such a manner that q = Q corresponds to the linear problem and q = 1 corresponds to 
the full quartic set of equations. Solutions of the full equation are thus smoothly connected 
to solutions of the simpler linear problem. 

There is no a priori way of knowing if all solutions have been determined, as well as 
which of these solutions are physically relevant. In the TICI2 formulation there are as 
many solutions as the number of equations; the lowest one corresponds to the ground 
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state of the system and the rest could be interpreted as monopole excitations or breathing 
modes. 

We have solved the non-linear system of equations for N =4, 16 and 40 bosons using 
the S3 interaction. In Fig. 1 are displayed the energies of two selected TICC2 solutions 
(only one for A'^ =4) and their evolution with the sequential scheme. These two solutions 
are represented by the solid and the long-dashed lines. The linear TICK solution is also 
displayed as the short-dashed line. Several comments are in order regarding these three 
solutions. 

(i) All of them tend to stability with increasing values of iVmax) and at iVmax = 30 they 
have already converged. 

(ii) There is a connection between the solid line and the TICK line in the quenching 
scheme and in both directions. This means that starting from TICK solution (q = 0) and 
slowly increasing the value of q one arrives to the solid line TICC2 solution when ^ = 1, 
and conversely, starting from the solid line TICC2 solution (g = 1) one arrives to the 
TICK solution by slowly decreasing q down to g' = 0. 

(iii) There is not a connection path which starts from the long-dashed TICC2 solution. 
When q is slowly decreased the solution disappears. 

In conclusion, even if the TICC2 solution represented by the long-dashed line has the 
lowest energy value, it should be considered as a spurious solution of the TICC2 method 
which has no relation with the physical problem. This discussion has been quite simplified 
because we have focused only on two cases. Actually, by generating at random starting 
points for the Newton-Raphson method, there appears a very large number of solutions 
which must be discriminated. 

The arguments of correctness of solutions, both to increasing N^nax or along the 
quenched path, do not really ascertain the goodness of a given solution. A precise but 
costly way is to compute the expectation value of the hamiltonian for any of the found 
solutions, but this is a hopeless task in the Fock representation when a large number of 
particles is involved. For this calculation one requires the matrix element 

(H) = ^ ^^ , (32) 



with a number of contractions scaling as N"^. 

An alternative way is to go to coordinate representation, where 



i<j ' i<j k<l 

+ ^EE'E%/«/™" + ---}'^Ho(ri,...,r^) (33) 

i<j k<l m<n 



In this expression, the correlation function / is defined as 



fir) = £ Sin) 2 [-^^] ''' L]l\\a\\ (34) 



n=l 
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and the sums over single pairs, double pairs, . . ., reflect the combinatorics of the selec- 
tion of annihilation operators for the contraction of the creation operators contained in 
: exp(S'^'^'*) :. The primes mean that no repeated particle indices appear in the multiple 
sums. The dots in Eq. (^) indicate up to N/2 sums. 

In Fig. 2 are displayed the pair correlation functions /(r) for the solutions considered 



in Fig. 1. With the corresponding wave function given in Eq. (33) we have computed 
the expectation value of the Hamiltonian by means of the Monte Carlo method with 
a Metropolis sampling, and the reader should be aware that the sampling with a wave 



function like (33) is very costly in terms of computation time. 

The simplest (and test) case corresponds to iV = 4. There we have not found spurious 
solutions, and the variational Monte Carlo (VMC) estimate for the energy corresponding 
to the lowest-energy TICC2 solution is —25.55 it 0.05 MeV, in a very good agreement with 
the normal TICC2 value, which is displayed in Table 1. 

The test has also been carried out for the normal and the spurious solutions of A^ = 16. 
For the solid line solution the VMC value is —1207 it 16 MeV, in good agreement with 
the normal TICC2 value, which is displayed in Table 1. However, the VMC value for the 
long-dashed solution turns out to be —661 it 171 MeV, with a very large variance, and 
very different from the corresponding TICC2 value of —1638 MeV. This is, in our mind, 
a definitive argument to consider the long-dashed line of Figs. 1-2 as a spurious solution, 
which should be ruled out. It is very interesting the agreement between TICC2 and 
TICC2-VMC calculations, and it is particularly stimulating if we realize that the TICC2 
method is not variational. The closeness of both results means that TICC2 determines 
appropriately the energy. 

The non-variational character of the SUB(n) truncations of CCM should not be consid- 
ered as a drawback. Variational methods are very popular in many-body calculations, but 
the only advantage of using them, apart from obtaining an upper bound, is just to secure 
an improvement of the ground-state energy when the Hilbert space is enlarged. On the 
other hand, the hierarchical structure of the SUB(n) truncation has some similarities with 
perturbation theory, in the sense that the exact value is not monotonically approached. 
Both methods may be quite far away from the exact (and unknown) ground-state energy: 
in Jastrow case because of using an inadequate correlation, and in CCM because of a too 
early truncation. 

It is worth mentioning that the problem of spurious solutions of the CCM method is one 



of the open questions related to this theory. In some simple cases |10], it has been possible 
to study the full map of solutions, but nevertheless these studies do not help too much 
in understanding the real origin of the spurious solutions. Here, the claimed advantage of 
the CCM, i.e. its non-linear character, turns out to be also a serious drawback. 

Once the TICC2 solution has been properly identified, it is time to comment on the 
physical results, going back to Table 1. To this end, we should have in mind that a 
nuclear bosonic system with the central S3 interaction tends to collapse when the number 
of particles A^ increases, so that it magnifies the role of the correlations. In Table 1 we 
can see that going from TICI2 to TICC2 results in an increasing gain of energy, ranging 
from 0.18 MeV for A^ = 4 to 963.24 MeV for A^ = 40. The importance of the non-linear 
terms contained in TICC2 increases with the number of particles. 
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It is also interesting to compare with the results provided by an hybrid ansatz which 
includes pair correlations as well as hyperradial excitations in a TICI2 scheme |Il|, 12 1. The 



combination of two-particle plus hyperradial correlations is the base of the hyperspherical 
harmonics description largely used in few-body systems |13, 14]. The entry TICI2+HR 



in Table 1 shows the results obtained by minimizing the Hamiltonian with respect to the 
correlation function /(/>, r), where p is the hyperradius. Such an ansatz is an obvious 
generalization of the TICK one. The comparison shows that the gain in energy due to 
hyperradial excitations is relatively small. TICI2 results change in less than 2% when 
hyperradial correlations are included. The gain due to TICC2 is much more important. 

Finally, to complete this analysis it is important to compare our results with those 
obtained using a Jastrow correlation function. The entry J-VMC in Table 1 shows the 
ground-state energies determined by a variational Monte Carlo procedure. For A^ = 16 and 
A^ = 40 a rough determination of the ground-state energy was made ||ll[] using a simple 
Jastrow correlation function fj{r) = 1 -|- oexp(— 6r^). For A^ = 4 we quote the Jastrow 
result given in Ref. Q. We can see that our TICC2 results are still much less bound 
than the Jastrow variational lower bound. This reveals the importance of three-body and 
higher order correlations. 

5 Final comments 

In this paper we have presented a way of determining the ground-state energy and wave 
function of a system of bosons described in terms of the translationally invariant formula- 
tion of the coupled-cluster method. The analysis has been limited to the case of two-body 
correlations (TICC2 approach), and has been carried out by using Fock space techniques. 
The correlation amplitudes satisfy a non-linear coupled set of equations, which has been 
explicated here, and has been solved for light and medium bosonic nuclear systems, con- 
sidering the S3 semi-realistic interaction of Afnan and Tang. 

To conclude we would like to mention the required further steps and the pending 
questions to be solved on line. The first obvious extension of TICC2 is to consider SUB (3) 
and presumably SUB(4) truncation, as it is suggested by the comparison with Jastrow 
results. Here the main problem to face up is to find a unique way of describing the 
translationally invariant cluster operators, both in coordinate and in Fock spaces. One 
may propose that the three-body operators should be a mixture of Ip-lh, 2p-2h and 3p- 
3h excitations, and the main questions are how to find a unique description (i.e. the object 
equivalent to S{n) of SUB(2)), and its relation to the characterizing quantum numbers. 

In previous studies restricted to light systems (or to heavier nuclei in the TICK approx- 
imation) the advantages of the coordinate representation have been stressed with regard 
to the actual computations and with respect to the richness of the correlation operators 
(i.e. the easy way to introduce spin and isospin dependence when dealing with true nu- 
clei). The form of the wave function is known [see Eq. (|33|)], and it seems simpler to 
determine a single function /(r) instead of a large number of numerical amplitudes S{n). 
The still open problem is to transform the Schrodinger equation with the ansatz (|30|) into 
an appropriate integro-differential equation for the two-body correlation function. 

The coordinate representation of the SUB(2) truncation of CCM is quite similar to the 
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familiar Jastrow form with only two-body correlations, where the model state is correlated 
by the action of the operator JIkj [^ + /(^«i)] • The expansion of this product generates all 
terms appearing in Eq. (^) and also additional terms which go beyond the independent- 
pair character of that equation. It is worth noting that a form like Eq. (|3^) was used time 



ago |15| to obtain hypernetted-chain equations for the case of state-dependent correlations. 
Moreover, the main advantage of the Jastrow ansatz with respect to the SUB(2) trunca- 
tion, particularly in the case of strong interactions, is its ability to tame the potential for 
all nucleon pairs. It has been suggested [^ to use this formal analogy to generate other 
truncation schemes for CCM which incorporate a Jastrow-like behaviour, and specially 
adapted to strongly short-distance repulsive interactions. 

The last open question is to find a TICC2 formulation for fermionic systems. Here 
again we expect serious algebraic difficulties related to the structure of the reference state, 
with several occupied shells and the corresponding complexity of the S^^'^^ operator. In 
our opinion, the formulation in coordinate space will simplify the problems related to the 
fermionic statistics. 
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Table 1: Binding energies (in MeV) of N bosons using the S3 interaction. 





iV = 4 


iV = 16 


A^ = 40 


TICC2 


-25.49 


-1234.86 


-8456.55 


TICI2 


-25.31 


-1130.94 


-7493.31 


TICI2 + HR 


-25.99 


-1145.11 


-7555.12 


J-VMC 


-27.15 


-1403±1.5 


-9570±7 



Figure captions 

Fig. 1. Energies (in MeV) of two selected TICC2 solutions (solid and long-dashed lines), 
plotted as a function of the number A'max of equations considered in the system dsH) . The 
TICK solution is also displayed as the short-dashed line. 
Fig. 2. Correlation functions /(r) as defined in Eq. (p^) for the solutions of Fig. 1. 
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